When is a bath a bath? Relaxation dynamics and thermahzation in a fermionic chain 
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We study thermalization in a one-dimensional quantum system consisting of a non-interacting 
fermionic chain with each site of the chain coupled to an additional bath site. Using a density 
matrix renormalization group algorithm we investigate the time evolution of observables in the 
chain after a quantum quench. For low densities we show that the intermediate time dynamics can 
be quantitatively described by a system of coupled equations of motion. For higher densities our 
numerical results show a prethermalization for local observables at intermediate times and a full 
thermalization to the grand canonical ensemble at long times. For the case of a weak bath-chain 
coupling we find, in particular, a Fermi momentum distribution in the chain in equilibrium in spite 
of the seemingly oversimplified bath in our model. 



Introduction. The time evolution of classical and 
quantum systems is deterministic. If a system in the 
thermodynamic limit reaches thermal equilibrium at long 
times, then we expect, however, that its physical prop- 
erties are determined by only a few parameters such 
as temperature, chemical potential, and pressure. This 
thermalization process is often studied in two different 
settings: (a) The system is in contact with a thermal 
bath, i.e., a large reservoir of thermal energy. The key 
assumptions commonly used in this setting are a weak 
coupling between bath and system and Markovian dy- 
namics, i.e., a very short correlation time in the bath. 
In this case the microscopic details of the bath become 
unimportant [TH5] ; for a simple example of classical ther- 
malization, see U]. (b) The system is closed, with parti- 
cles being able to exchange energy and momentum among 
each other, so that the closed system can explore phase 
space, constrained only by conservation laws such as to- 
tal energy and particle number. An important difference 
between the two scenarios is that in the first case tem- 
perature, chemical potential and pressure are parameters 
determined externally by the bath. In the latter case, on 
the other hand, these parameters are Lagrange multipli- 
ers fixing the values of the conserved quantities [5, 6J. 

In this letter we want to study these two settings simul- 
taneously using a model which can be either viewed as a 
closed quantum system or as a chain coupled to a simple 
bath. Thermalization, in both cases, requires: (I) observ- 
ables become time independent and all currents vanish 
{equilibration); and (II) time averages can be replaced 
by statistical averages over ensembles with a restricted 
number of intensive parameters [35] , and are indepen- 
dent from initial conditions (ergodicity) [?]• The rather 
old but fundamental problem of non-equilibrium dynam- 
ics and thermalization in closed quantum systems has 
been put again into focus by experiments on cold quan- 
tum gases which are very well isolated from their sur- 



roundings jHHHli as well as by the development of new 
numerical techniques to study dynamics in many-body 
systems [T2l - [l9] . This has lead to numerous simulations 
of non-equilibrium dynamics in closed quantum models 
where the question whether or not thermalization occurs 
has not always been easy to answer due to the finite nu- 
merical simulation time |20H23| . 

Closed quantum systems. The time evolution of an 
initial state |^o) = — 0)) is unitary and given by the 
Schrodinger equation. Therefore |^'(t)) remains a pure 
state for all times t. Since ensemble averages describe 
mixed states such a description cannot apply to a finite 
closed quantum system as a whole. Only a subsystem can 
be in or close to a thermal state with the rest of the sys- 
tem acting as an effective bath. Furthermore, contrary 
to a classical system, every quantum system has expo- 
nentially many conserved quantities, e.g. the projection 
operators P„ = |£'„) (i?„ | onto the eigenstates of a system 
with discrete spectrum, \E„) ^Bj. However, it is usually 
assumed that only the local conserved quantities are of 
relevance for thermalization. A local conserved quantity 
can be represented for a lattice systems as Q,n — J2j 
where is a density operator acting on lattice sites 
J, + 1, • • • ,j + m with m finite. Here we want to con- 
centrate on the case of generic one-dimensional quantum 
systems with a small number of local conservation laws, 
i.e. the total energy and particle number. Thermaliza- 
tion in closed integrable models, where the number of lo- 
cal conservation laws increases linearly with system size 
\24\ [25] , has been investigated with the help of numerical 
simulations in recent times as well [26.-29 . 

We consider the non-equilibrium dynamics ensuing af- 
ter preparing the system in a pure state |^o) which is 
not an eigenstate of the Hamiltonian. Using a Lehmann 
representation we can write |^o) = X^n^"!^") where 
\En) are the eigenstates of the Hamiltonian H the sys- 
tem evolves under. Furthermore, we restrict ourselves to 
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typical states with a macroscopic number c„ ^ |37j . We 
can now easily calculate the long-time mean 

O = lim V- /'dte'(^'"-^")*<c„(i?„|0|i?„) 

n 

of an observable O, where we have set 1. The second 
line of ([T|) is often called the diagonal ensemble. Here we 
have assumed that the system is generic, i.e., that degen- 
eracies play no role. If the observable becomes stationary 
at long times then its value Ooo — ^^i^t'^oo{ilj{t)\0\ip{t)) 
has to be equal to the long-time mean, Ooo = O. Note 
that this is only possible in the thermodynamic limit. 
Otherwise observables show revivals on timescales of the 
order of the system size. Taking the thermodynamic limit 
is thus essential; a finite system can never thermalize. 

If a subsystem of an infinite system containing the 
observable O equilibrates and the value Ooo does not 
depend on details of the initial state, then the remain- 
ing open question is which ensemble describes the equi- 
librated system. If we have two statistically indepen- 
dent subsystems A and B then the density matrix p of 
the whole system is given by p — pA® Pb, thus Inp = 
In 0ln pb- Second, the density matrix itself should be- 
come time independent once the system has equilibrated 
and the von-Neumann equation implies p = —i[H, p] = 0. 
Thus the general density matrix under consideration has 
to be of the form p = exp(— A„Q„)/Z where Q„ are 
the conserved quantities of the system [S]. The partition 
function Z is a normalization factor such that Tr p = 1 . 
We stress again that the intensive parameters A„ are not 
given externally but rather are Lagrange multipliers de- 
termined by the set of equations 

(*o|Qn|*o> =Tr{Q„p}. (2) 

If we include all projection operators into our density 
matrix, Q„ = P„, then it follows immediately from ([2| 
that (0)p = Tr{pO} is identical to the diagonal ensemble 
as given in Eq. ([T]) Having to use infinitely many 
Lagrange multipliers is expected because is always 

a pure state and the system as a whole therefore does not 
thermalize, because it does not fulfill condition (II). 

In this letter we focus on the generic situation where we 
split our system S = AuB into a bath B and a subsystem 
A and consider observables acting only on subsystem A. 
We concentrate on the following questions: How does a 
subsystem A without intrinsic relaxation processes equi- 
librate when coupled to a strongly correlated but simple 
and possibly non-Markovian bath B7 Which statistical 
ensemble gives the expectation values of observables in 
A in the equilibrated state? 

Model Hamiltonian. To investigate some aspects of 
the questions raised above we consider a simple model 
system with Hamiltonian |30| 
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E {'h - 1/2) - 1/2) . (3) 

The first term describes a chain of free spinless fermions 
with hopping amplitude J and is the subsystem A we 
study the thermalization of. The 'bath' B consists of 
extra sites, coupled to the chain sites via a hybridization 
7 (second term), and we also include a nearest-neighbor 
interaction Vs between the bath sites (third term). 

As initial states for the time-evolution with the Hamil- 
tonian ([3| we will consider, on the one hand, the 
ground state 1^-^,(70,70)) = |^'(Jo,7o,K = 0))o of 
the non-interacting model with hopping parameters Jq 
and 7o as well as the ground state |5'o^( Jo, Jq, 70)) = 
|\l/( Jo, Jq, 7oi = 0))o of ([3]) with an additional hop- 
ping Jq between the bath sites. In order to study 
the time evolution under the interacting Hamiltonian 
we use a time-dependent density-matrix renormalization 
group (DMRG) algorithm [31], a method which has al- 
ready been applied to study other one-dimensional mod- 
els [21^,22, 32 . We choose open boundary conditions with 
a chain length of L = 51. The number of states kept 
in the truncated adaptive Hilbert space varies between 
X — 400 and x — 800. For a global quench as considered 
here it is well known that the entanglement entropy be- 
tween two subsystems increases linearly with time. Since 
the maximal entanglement which can be represented in a 
truncated Hilbert space is limited by lux, there is a max- 
imum time iinax up to which we can reliably simulate the 
time evolution. For the cases considered here this time 
scale is given by Ji,nax ~ 15 — 25. 

Results. First, we will concentrate on the relaxation 
dynamics at low particle densities. As an example, we 
show in Fig. [1] results for a quantum quench with = 11 
particles. Shown are results for the one-point correlation 
functions 

C,{t) ^ {C,)t = (*(t)|c;^^,)/2C(L+i)/2+,l*W> • (4) 

In all correlation functions oscillations with a character- 
istic frequency are visible. These oscillations can be un- 
derstood from an equation of motion approach. We de- 
fine the three time-dependent expectation values fq{t) = 
{^{t)\c\c,\^{t)), g,{t) - {^{t)\s\s,\^>{t)) and p,{t) = 

(vl'(t)|4,s,|*(<)) where c, = V2/(i + 1) sin(gj)c, 
with allowed momenta q — n-n /{L + 1), n = 1, . . . , L. 
Then, using Heisenberg's equation of motion, a Hartree- 
Fock decoupling of the quartic terms, and the additional 
assumption of an instantaneous dephasing |33| . we find 
the following system of coupled equations 

fqit) = -9q{t) = 27^g(i) , = Bg{t)R,{t) , 

R,it) = -iim - gM - B,{t)r,{t) , (5) 
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FIG. 1: (Color online) Cj{t), Eq. Q, for j = 0, . . . , 5 (top to 
bottom) for a quench with initial state |^'o(l, 1)), and Hamil- 
tonian H with J = 1, 7 = 1, and Vs = 1 a.t low densi- 
ties. DMRG results (symbols) are compared to the solution 
of Eq. m (lines). 



rq{t) = Repq{t), Rq{t) = Imp,(i), 

Ea + 2Vs/{L + I)[c0s2((?)5,_,(t) - 



with Eq = — 2cosg 
and Bq{t) = -14 - 

sin^(g)5q(t) + X]fc(l-cosfccos(j)5fe(t)] « -Vg-Eq = B, 

We solve the set of equations (|5| numerically and the 
results up to intermediate times are in excellent agree- 
ment with the DMRG data, see Fig. [T] Using further ap- 
proximations, we analytically find that the oscillation fre- 
quency is given by fl^ = Bl + {2^f with f]2_^o ss l + {2'yY 
and depends only weakly on q [33]. This means that the 
dephasing process is very slow. For longer times and 
short distances we see that the amplitude of the oscil- 
lations in the DMRG data is decaying faster than pre- 
dicted by our equations of motion approach. Here it is 
important to realize that due to the Hartree-Fock de- 
coupling the equations of motion effectively describe the 
time evolution under a free particle Hamiltonian. This 
approach therefore takes only the slow dephasing process 
discussed above into account. The additional decay seen 
in the DMRG data is due to slow relaxation processes 
involving energy-momentum transfer between interact- 
ing particles which are not captured in our equations of 
motion approach. 

A much faster relaxation occurs if we increase the par- 
ticle density with a maximum in the relaxation rate at 
half-filling. The DMRG data for a quench in the half- 
filled case in Fig. |2]show indeed that the system almost 
completely equilibrates within the simulation time tmax- 
Due to particle-hole symmetry of the Hamiltonian and 
the initial state we have Co{t) = 1/2 and C2j{t) = 0. For 
odd distances we now see, instead of long-time oscilla- 
tions, an exponential damping which allows us to extrap- 
olate the correlation functions and to read off the value 
for Cj {t 00) L33j . Due to the lightcone-like spreading of 
correlations [T71I3^ . the short-range correlation functions 
in the middle of the chain are, for the time range shown 
in Fig. |2] not affected by the boundaries and are almost 
indistinguishable from those for an infinite system. By 
extrapolating our numerical data we thus approximately 
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FIG. 2: (Color online) DMRG data (symbols) for Cj(t) at 
half-filling for a quench with initial state l^fo' (1,0.6, 1)), and 
Hamiltonian H with J = 1, 7 = 1, and Vs = 1- The lines are 
the thermal expectation values {Cj)T- 
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FIG. 3: (Color online) (a) fg{t = 0) for |*;j^(l, 0.6, 1)) (cir- 
cles), fq{t — 5) (triangles) and Fermi function fit {/q)T=o.7J, 
and the extrapolated distribution fq{t — >■ 00) (diamonds) 
with a fit {fq)T=j. (b) fq{t — > 00) (diamonds) compared 
to the thermal average TT{fqe-"^'^}/Z (solid line) and {/,)t 
(dashed hne) where T/J = 0.54 is fixed by Eq. (c) fq{t) 
for q — 137r/52 (solid line) and q = 167r/52 (dashed line). 



obtain C'j{t — >■ 00) for a system in the thermodynamic 
limit. 

The corresponding distribution function fq{t), shown 
in Fig. [sja), has already become completely smooth af- 
ter a short time, Jt — 5, and can be well fitted by a 
free fermion distribution function {fq)T = l/{e^''^'^ + 1). 
However, the system has not fully equilibrated yet. 
Fig. ^c) shows that we have two distinct relaxation 
regimes. In regime Rj we have a relatively quick reshuf- 
fling in the distribution leading to a prethermalized state 
[HI |3S] . This is followed by a slow drift of the occupation 
numbers in regime Rjj which, when extrapolated in time. 
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FIG. 4: (Color online) The initial distribution (circles), 
fq{t — )■ oo) from DMRG (diamonds), the thermal distribution 
Tr{fqe~^''^}/Z (solid line), and the free fermion distribution 
{fg)T (dashed line) with: (a) 7 = 1, T/J = 0.19, (b) 7 = 0.6, 
T/J = 0.18, and (c) 7 = 0.2, T/J = 0.33. 



leads to the final distribution for the equilibrated state. 
While both distributions can be well fitted by {fq)T, the 
temperature should not be used £lS £L fitting parameter 
but rather be determined by energy conservation. We 
therefore expect that the equilibrated system is described 
by the ensemble, p = exp(— i7/T)/Z, with chemical po- 
tential fi — due to particle-hole symmetry. The tem- 
perature is determined by Eq. (|2| with Q„ replaced by 
H. The Ihs of (|2| is now an expectation value for a non- 
interacting system and can be obtained analytically. The 
thermal average on the rhs is calculated using a static 
DMRG calculation [31 with x = 50 states kept. For the 
particular quench in Fig. [s] we find T/J = 0.54. This 
then allows the calculation of {Cj)T = Tr{(7je^^/-^}/Z 
by DMRG shown in Fig. |2] The results for the corre- 
sponding distribution function are shown as a solid line 
m Fig. ^h) and agree well with the time extrapolated 
values, demonstrating a local thermalization. 

If the additional sites are to represent an effective 
bath, the distribution function in the chain should be- 
come a Fermi distribution. However, as can be seen in 
Fig. [sj^b), (/g)T=o.54J differs significantly from the equi- 
librium distribution. One obvious reason is that the effec- 
tive coupling between chain and bath in the thermal state 
~ 7(s|ci)T=o.54J ~ 0.287 is not small. Next, we there- 
fore consider cases where we successively reduce the cou- 
pling 7. In order to be able to still find the equilibrated 
state within the limited simulation time we now use as 
initial state |5'q(1,7)) which yields a much smoother ini- 
tial distribution. Results for different coupling strengths 
7 are shown in Fig. |4] We indeed find that the mo- 
mentum distribution in equilibrium now approaches the 
free fermion distribution with the temperature deter- 



mined by Eq. ([2|. At 7 = 0.2 the effective coupling 
between chain and bath cx 7(sJci)T=o.33J ~ O.O67 is very 
small. Apart from the usual Pauli blocking there is an- 
other mechanism which explains the very weak coupling 
between the subsystems. Because the nearest neigh- 
bor occupation (7t.^7t.^j^)t=o.33J — 0.1 is also small we 
can approximately project out all states where nearest- 
neighbor sites in the bath are occupied. This leads to 
an effective density-density interaction c>c (7^/14)71^^71^ 
between the subsystems, leading to a slow relaxation 
for small 7 In this strong coupling limit, the hy- 

bridization part of the Hamiltonian ([3| also gets pro- 
jected cx '^{s^jCj + h.c.){l ~ n^_i){l — n^_^_i) explaining the 
small value for the effective coupling given above. Thus 
the interactions help to decouple the two subsystems ex- 
plaining the almost perfect free Fermi distribution in the 
chain for 7 = 0.2. 

Conclusions. We have studied thermalization in a 
strongly correlated model which can be viewed either as 
a closed quantum system or as a free fermionic chain 
coupled to a bath. Contrary to the common approach 
of using a Lindblad equation to study open quantum 
systems, our model has a microscopically specified bath. 
Therefore we can simulate the non-equilibrium dynamics 
of system and bath and directly compare the two differ- 
ent viewpoints. For low particle densities we have shown 
that an equation of motion approach on the Hartree-Fock 
level is sufficient to quantitatively describe the interme- 
diate time dynamics. At this level only slow dephasing 
processes are captured. For the future it seems promis- 
ing to use a higher order decoupling which might also 
capture the faster relaxation processes which we observe 
in the numerical simulations. While the relaxation rate 
F ~ Vs{nfnf_^_i) at small interactions or low densities is 
too small to observe equilibration within the limited nu- 
merical simulation time we do observe thermalization at 
stronger interactions near half-filling where F is larger. 
We note that the relaxation rate changes continuously 
with the microscopic parameters of the model so that the 
definition of a 'non-equilibrium phase transition' based 
on the accessible simulation time tmax seems problematic 
|21| . Most interestingly, we find that strong interactions 
lead to an effective disentanglement between the subsys- 
tems and therefore increase decoherence times. Further- 
more, even an extremely simple bath where Markovian 
dynamics can not be taken for granted can be sufficient to 
fully equilibrate a subsystem without intrinsic relaxation 
processes. 
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the National Natural Science Foundation of China (Grant 
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Equation of motion approach 



Using Heisenberg's equation of motion. 



We set up a set of equations for the following three 
time-dependent expectation values 



Ut) - (vl/(t)|ctc,|*(t)) , 
g,it) = (vl/(t)|4s,|M/(t)) , and 
p,it) - {^it)\cls,\^{t)) . 



(6) 



6 = i[H, O] 



(7) 



with the Hamiltonian given by Eq. (1) of the main text 
and the standard fermionic commutation relations, one 
finds, with — s^jSj, 



m 

9S) 



2-fR,{t) , (8) 



L-1 L 

^ Sin [qj] {^{t)\cls,nf+,\^{t)) + ^ sin [qj] {^{t)\clsjnf,,\^{t)) 



L + 1 



J=2 



We have introduced the real functions (t) = Re pg {t) 
and i?q(t) = lmpq{t). This set of coupled equations is 
exact. 

To solve this set of equations we apply two approxima- 
tions. Firstly a Hartree-Fock decoupling is used, i.e. we 
apply Wick's theorem so that we have only two point 
correlation functions present, for example 

One should note that, amongst other effects, this approx- 
imation leaves 



(10) 



and hence gq{t) + fq{t) = Nq becomes independent of 
time. Therefore, by performing the Hartree-Fock decou- 
pling, we lose all relaxation processes which can reshuffle 
the occupation of the momenta. Nonetheless, as demon- 
strated in Fig. 1 of the main text, for small densities 
this approximation is sufficient to get good quantita- 
tive agreement with DMRG calculations for intermedi- 
ate times. The second approximation is "instantaneous 
dephasing", which means that all off-diagonal elements 
of {'i' {t)\c^qCk\'i {t)) , etc., are taken to 'instantaneously 
dephase' and we keep only diagonal terms. For a peri- 
odic system this would be guaranteed by translational 
invariance, here it amounts to disregarding finite size ef- 
fects from the boundaries. This means that all two point 
correlation functions are diagonal in momentum space. 
Following this we have 

Mt) = -9qit) = 27-Rg(i) , rq{t) = Bq{t)Rq{t) , 
R,{t) = -iW) - gM ~ Bq{t)rq{t) , (11) 
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FIG. 5: Main: A quench with initial state |*I'o(l,l)), and 
Hamiltonian H with J = 1, 7 = 1, and Vs = 1 at low densities. 
Shown is fq{t = 0) (circles), fq within Hartree-Fock (squares), 
and fq obtained by DMRG (diamonds). Inset: Q.q obtained 
by DMRG (symbols) and within Hartree-Fock (line). 



with 



2Vs 
L+1 



cos^ [q] g^-q{t) ~ sin^ [q] gq{t) (12) 



- ^ (1 - COS [k] cos [q]) gk{t) 



The coupled first order differential equations, given by 



Eq. (|ll|), can then be solved iteratively. The first line in 
Eq. 



12) is a 0{1/L) finite size correction. 



The oscillation frequency of Rq{t) is the same as that of 
fq{t) and gq{t) and can be extracted from these equations 



7 



analytically. We can write a second order differential 
equation for Rq{t): 

R,{t) + {A^^ + Bl{t)) R,{t) = -B,{t)r,{t). (13) 



tion. such that Bq ^ 



One finds, with a weakly time dependent bath occupa- 
Rg{t) + {Aj' + B'^it)) Rq{t) = 0. (14) 



For small particle densities in the bath we can approxi- 
mate Bq{t) ^ —Vs — Eq = Bq which glvcs 



4^2 



(15) 



The Hartree-Fock decoupled solutions oscillate with the 
frequency f2g which is only weakly g-dependent, see inset 
of Fig. [5] This explains why no dephasing effects are seen 
in the Hartree-Fock solution on the timescales that we 
consider, see Fig. 1 in the main text. The approximation 
Bq{t) ss Bq is ensured in our case by the small bath 
occupation. For example, in the initial state |\1/q(1,1)), 
at a density of 0.11 particles per site, we have 



L 



0) ^ 0.0344< (16) 



and 



L 



0)cosg = 0.0314< (17) 



Contrary to the Hartree-Fock results, the DMRG data 
show an additional relaxation, see Fig. 1 of the main text. 
A signature of the beginning of this relaxation can also be 
seen in the long time mean of the distribution function, 
/g, see Fig. [s] which shows a redistribution of the occu- 
pation of quasi-momenta around the Fermi momentum. 
The low density relaxation rate F ^ VsJ2j{''^f''^f+i)/-^' 
however, is small so that we can not see full thermaliza- 
tion within the DMRG simulation time imax- 
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FIG. 6: Momentum space extrapolation for a quench with 
initial state |4'o(l,7)), and Hamiltonian H with J = 1 and 
Va = 1 where (a) 7 = 1, and (b) 7 = 0.2. Shown are the mo- 
menta q = lln/{L + 1) (upper curves) and q = 21tt/{L + 1) 
(lower curves). The dynamics (symbols) are compared with 
the fit (dashed line), and the thermal average at the appro- 
priate temperature (solid line). Fitting is performed for times 
greater than the solid vertical line. The arrows on the right 
hand side show fq{t — >• 00) from (a) Eq. (211, and (b) Eq. (22 1. 
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FIG. 7: Real space extrapolation for a quench with initial 
state |^o'(l, 0.6, 1)), and Hamiltonian H with J = 1, 7 = 
1, and Vs = 1. The DMRG data (symbols) are compared 
with the fit (dashed line), and the thermal average at the 
appropriate temperature (solid line). Plotted is Cjii) with 
(a) J = 1, and (b) j = 3. Fitting is performed for times 
greater than the solid vertical line. The arrows on the right 
hand side show Cj(t — > 00) from Eq. (|23[). 



Particle— hole symmetry 

We define iJ' as the Hamiltonian given by Eq. (1) in 
the main text, but with hopping in the bath included: 



L-l 



(18) 



i/', and therefore also all half-filled groundstates, has 
particle-hole symmetry. The Hamiltonian B' is invariant 
under the mapping 



^ (-l)^'+is] 

^ (-1)^"+^., 



(19) 



which exactly describes particle-hole inversion. 

In our analysis we have considered two correlation 
functions. Firstly the real space two point correlation 
function Cj{t), defined by 



(20) 
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FIG. 8: Distribution function fg in the chain for a system 
of size L = 51. = 0) is shown for the initial state 

|*"(1,0.6, 1)) (circles), and the initial state |*o(4.48, 0.8)) 
(squares). In both cases we time evolve with the same Hamil- 
tonian H with J = 1, 7 = 1, and Vs — 1. fq{t — >■ 00) after 
time evolving |^'o^) (diamonds) and l^&o^) (triangles) are com- 
pared with the thermal average TT{fqe~^^^} /Z (solid line) 
where T/J = 0.54 is fixed by the initial energy, see text. 
Alternate q points are plotted for the two quenches to aid 
clarity. 



(191 



H - 

:l/2, 



' H 
and 



for time evolution with H given by Eq. (1) in the main 
text. Under the mapping given by Eq 
and Y^q) — )■ Iv&o) and one finds that Cj=Q{t) 
C2j{t) = for non-zero j. Secondly we analyzed the 
momentum distribution in the chain, fq{t). For the initial 
states under consideration, and therefore for all times, 
this can be shown to satisfy fq{t) + /^_g(t) = 1 by using 
the same mapping. 



niq = ll7r/(i + l)) « 2.7 and n{q = 2l7r/(L + l)) « 1.9. 
The relaxation rate, r{q = Il7r/(L + 1)) « 6.7 • lO'^ and 
T{q = 2l7r/(L + 1)) w 2.7 • 10"^^ is of the same order of 
magnitude for all momenta hinting at one dominant re- 
laxation process. Let us stress that in this case the result 
for fq{t — > 00) depends only weakly on the extrapolation 
procedure used, i.e. time averaging or fitting with differ- 
ent fit intervals, with a variation in fq{t — >■ 00) which is 
about the symbol size used in the corresponding plots in 
the main text. 

For 7 — 0.2, see Fig. |6jb), such a simple fitting func- 
tion will no longer work due to the presence of various 
oscillation frequencies. A Fourier analysis confirms that 
there is more than one oscillation frequency involved, but 
the times are not sufficient to extract how many there are 
and what their magnitudes may be. Instead we trace out 
the overall trend by fitting to 



/,(Ji»l)~/,(t^CX))+, 



-rt 



(22) 



F captures the gradual drift of the oscillations which can 
also be seen by using running averages. This procedure is 
robust when choosing a variety of different time regions 
over which to perform the fitting, and gives again errors 
smaller than the symbol sizes used in the plots of the 
main text. 

For quench II we extrapolate in real space. Fig.|7]shows 
the fitting for Ci (t) and C3 {t) . We fit the dynamics with 

Cj{Jt > 1) ~ Cj{t -^00)+ ae"" (cos [fit - <j)] + b) (23) 

As for quench I, only times right of the vertical lines in 
Fig. |7] are used for fitting. 



Fitting and extrapolation 

In this section we explain the fitting and extrapolation 
techniques we used to find our long time data at half- 
filling. For quench I we extrapolated in momentum space. 
The long time limit for fq and 7 = 1 and 7 = 0.6 can be 
found directly by time averaging the data or by fitting to 

fq{Jt > 1) ~ fq{t ^ 00) + ae"" cos [nt - (j)] . (21) 

This functional form takes into account exponential re- 
laxation and a simple oscillation, see Fig. ^a.). Note 
that the fitting is only performed right of the solid line at 
Jt « 5, to ignore the effect of the short time dynamics. A 
Fourier analysis confirms that there is one dominant fre- 
quency in the dynamics of fq(t). However, in contrast to 
the low density case this frequency is momentum depen- 
dent. In particular for the case 7 = 1 shown in Fig. l6ja). 



Initial state independence 

After relaxation the equilibrium state should depend 
only on the energy in the system. As example, we take 
two initial states, and |^o^), constructed to have the 
same energy after a quench 

E = {^'^\H\^'^) = (^l\H\^l) . (24) 

The two different initial states are time evolved with the 
same Hamiltonian, with J = 1, 7 = 1, and Vs = 1. We 
use the initial states |\E'q^(1, 0.6, 1)) (same as in Fig. 3 
of the main text) and |5'q(4.48, 0.8)) with energy E = 
—51.19. Fig. [8] demonstrates that both states evolve to- 
wards the same equilibrium state, well described by the 
grand canonical ensemble Tr{/ge~^/-^}/Z with the tem- 
perature T/J = 0.54 fixed by = Tr{He-"/'^}/Z and 
fi — due to particle-hole symmetry. 



